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Abstract. We propose a theory of ripples in suspended graphene sheets based on 
two-dimensional elasticity equations that are made discrete on the honeycomb lattice 
and then periodized. At each point carbon atoms are coupled to Ising spins whose 
values indicate the atoms local trend to move vertically off-plane. The Ising spins are 
in contact with a thermal bath and evolve according to Glauber dynamics. In the 
limit of slow spin flip compared to membrane vibrations, ripples with no preferred 
orientation appear as long-lived metastable states for any temperature. Numerical 
solutions confirm this picture. 
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1. Introduction 

The experimental discovery of graphene (single monolayers of graphite) [U |2] and its 
extraordinary properties have fostered an enormous literature reviewed by a number 
of authors [3[ HI [5l [6j. The first observations of suspended graphene sheets showed 
evidence of ripples: nanometric ondulations of the sheet with characteristic amplitudes 
and wave lengths that do not have a preferred direction [7J [8] . Ripples are expected to 
have a significant impact on electronic transport in graphene [9] . There have been some 
attempts to explain rippling as an instability or an equilibrium phase transition. For 
example, at zero temperature, electron-phonon coupling may drive the graphene sheet 
into a quantum critical point characterized by the vanishing of the bending rigidity of 
the membrane [10]. Ripples arise then due to spontaneous symmetry breaking in a 
buckling transition [10J. The connection between rippling and electronic properties was 
also investigated earlier [HI [12]. Monte Carlo calculations also showed ripple formation 
at different temperatures and evidence that ripples may be connected to variable length 
cr bonds of carbon atoms [13J. The possibilities that the ripples are due to thermal 
fluctuations [HI [14] or to adsorbed OH molecules on random sites [15] have also been 
explored. There is active research about the effects of ripples and strain on electronic 
properties, including possible strain engineering [IB] ; see the review in [6j. 

Recent experiments have shown both micron long ripples on graphene p2] and 
very short Angstrom-sized ripples in bilayer graphene [18]. On large graphene sheets 
suspended on substrate trenches Bao et al have shown that long ripples along the 
direction of applied stress can be thermally induced and controlled by thermal cycling 
in a clean atmosphere [17] . In contrast to [7], rectangular graphene sheets on substrate 
trenches are clamped only on two of their sides and free on the other two [IT]. The 
resulting ripples are much longer (wavelength about 5 microns) than those observed on 
samples pinned on all their sides [TJ E]- Bao et al interpret their experimental results by 
means of a theory of wrinkling in thin elastic sheets based on simplified stationary Foppl- 
von Karman equations [E]. A different question is which ripple solutions are selected 
by the dynamics in the sheet. Our numerical simulations of the time-dependent von 
Karman plate equations (discretized on the honeycomb lattice) show that a suspended 
sheet that is initially in a rippled state remains rippled whereas a flat sheet remains 
flat. In the absence of longitudinal stretching, the initial condition selects whether the 
graphene sheet remains flat or it acquires curvature or ripples. 

In bilayer graphene, Angstrom-sized ripples seem to be caused by the preparation 
of the sheet. Some solvent gets trapped between the two layers of the graphene bilayers 
and the stresses arising during drying provoke these extremely short ripples [18] . This 
view is supported by the fact that these ripples do not appear in similarly prepared 
monolayer graphene [H]. Angstrom-sized ripples can be described by a two-dimensional 
(2D) Ising model with antiferromagnetic couplings involving nearest and next-nearest 
neighbors [20]. The Ising spins of [20] are not coupled to the membrane displacement 
vector at each site. 
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In this paper, we explain the small-scale ripples observed in suspended monolayer 
graphene sheets IE] differently. Firstly, observation of ripples with a high-resolution 
transmission electron microscope (HRTEM) implies that the graphene sample is 
bombarded by a low-intensity electron beam. Even though the electron beam does not 
have sufficient energy to knock off carbon atoms, it may create defects by inducing 
bond rotation and certainly excite the atoms. Thus an observed graphene sample 
is continuously excited and cannot be considered to be in equilibrium. We suggest 
that the electron beam may push the carbon atoms vertically away from the planar 
configuration of the graphene sheet in a random fashion. As a possible model we consider 
a membrane with atoms coupled to Ising spins whose sign indicates the atoms trend to 
move upward or downward from the planar configuration. These Ising pseudo spins do 
not correspond to true spins but they exert forces at each atom site on the graphene 
sheet and interact with each other through their coupling to the membrane. Related 
Ising spin-membrane models with antiferromagnetic coupling have been proposed to 
explain wrinkling in membranes associated with partially polymerized vesicles [22J. 
Elasticity coupled differently to a soft spin has been used to study two-dimensional 
melting mediated by dislocations [23J. In a simpler context, we have found out that 
one-dimensional mechanical systems coupled to Ising spins that flip randomly according 
to Glauber dynamics [21] have a buckling equilibrium phase transition at a critical 
temperature that increases with sample size. More importantly, these systems exhibit 
ripples at any temperature in the limit as the spin relaxation time is much longer 
than the vibration periods of the mechanical system [21]. The ripples are metastable 
quasi-equilibrium states about which the mechanical system experiences small and rapid 
oscillations. 

In the case of suspended graphene sheets, we can describe the position of carbon 
atoms by discrete elasticity equations previously used to model defect dynamics in 
graphene [23 123 123 EE]- Neighboring carbon atoms in a graphene sheet attach to 
each other using three of their bonds. The fourth bond is not saturated and, similarly 
to the effect of the electron beam, may try to pull the atom upward or downward from 
the flat sheet configuration. As we indicated before, this trend is modeled by coupling 
each atom to an Ising spin. These spins are in contact with a thermal bath and flip 
randomly according to Glauber dynamics. After a transient stage, a random initial spin 
configuration gives rise to a number of spin-up and spin-down domains that, in turn, 
produce domains of atoms displaced upward or downward from the planar configuration. 
These randomly oriented domains form the ripples. Once they appear, the ripples change 
slowly through spin flips of the atoms in their boundaries followed by changes in their 
vertical displacement. 

The rest of the paper is as follows. The model we use is described in Section 
[21 In Section |3j we show that the model has a buckling equilibrium phase transition 
at a critical temperature 9 C which turns out to be many times higher than room 
temperature for the samples we consider. For 9 > 9 C , the spins are unpolarized 
and the planar graphene sheet is thermodynamically stable. Below 9 C , the spins 
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acquire a finite polarization and the graphene sheet buckles. This phase transition 
corresponds to a bifurcation of the macroscopic equations of motion for mean values of 
the elastic displacements and the spin polarization. Since TEM observation results in 
a low-intensity but continuous bombardment of the graphene sheet by electrons, true 
thermodynamic equilibrium is never reached. In Section HJ we discretize the equations 
of motion on the honeycomb lattice, write them in primitive coordinates and periodize 
them, so that integer multiples of displacements on primitive directions leave the lattice 
unchanged. We then solve numerically the resulting stochastic equations of motion in 
the limit of slow spin flipping and fast membrane vibrations. We show that long-lived 
ripples with no preferred orientation emerge after a short transient even from random 
initial conditions. These ripples correspond to those observed in experiments E] . We 
also study how the presence of defects in graphene affects rippling. Lastly, Section [5] 
contains our conclusions. 



In the continuum limit, elastic deformations of graphene sheets have the free energy of 
a membrane 



where (ui,u 2 ) = (u(x,y),v(x,y)), w(x,y), k, A and \i are the in-plane displacement 
vector, the membrane vertical deflection, the bending stiffness (measured in units of 
energy) and the 2D Lame coefficients of graphene, respectively [29]. V = (d x ,d y ) is 
the 2D gradient and V 2 the 2D laplacian. In (J2J) we have ignored the small in-plane 
nonlinear terms d x .ud x .u + d x .vd x ,v. 

Carbon atoms in the graphene sheet have a bond orbitals constructed from sp 2 
hybrid states oriented in the direction of the bond that accommodate three electrons 
per atom. The remaining electrons go to p states oriented perpendicularly to the sheet. 
These orbitals bind covalently with neighboring atoms and form a narrow tt band that 
is half-filled. The presence of bending and ripples in graphene modifies its electronic 
structure. [5] Out-of-plane convex or concave deformations of the sheet have in principle 
equal probability and transitions between these deformations are associated with the 
bending energy of the sheet. A simple way to model this situation to consider that 
out-of-plane deformations are described by the values of an Ising spin associated to each 
carbon atom. We shall assume that the spins interact with each other only through 
their coupling with the membrane and have an energy 



where / has units of force per unit area. At any time t, a spin located at (x, y) may 
change its sign according to Glauber's dynamics: the atom-spin system may experience 



2. Model 




(1) 



(2) 




(3) 
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a transition from (it, v, w, er) to (u, v, w, R^ x y ^cr 



W, 



(T\U, V, W 
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2 L 



at a rate given by [2T] 

P(x,y)a(x,y)] , 



j3(x, y) = tanh 



fa 2 w(x,y) 
9 



(4) 
(5) 



Rr Xj y^<T is the configuration obtained from cr = {a(x,y)} (for all points (x,y) on the 
hexagonal lattice) by flipping the spin at the lattice point (x,y). The parameter a 
gives the characteristic attempt rate for the transitions in the Ising system. Since the 
bending energy of the graphene sheet is k, the attempt rate should be proportional to 
an Arrhenius factor, a = aoe~ K ^ e , where ao is a constant. At room temperature, this 
is a large factor, 10 -15 , which bridges the gap between the phonon time scales and the 
observed damping time for lattice defects (a few seconds [30J ) . 

The total free energy F = F g + F s provides the equations of motion 

\Vw\~~ 



p 2 d t u = X d x \ d x u + d y v + 

+ fi d y (d y u + d x v + d x wd y w) 



p 2 d\v = X dy 



d x u + d y v + 



\Vw\ 



+ n d x (d y u + d x v + d x wd y w) 



p 2 d 2 t w = X d x 



+ X8 V 



d x u + d y v + 



d x u + d y v + 
Vwl 2 



\Vw\ 



dyW 



+ pd x [2d x u+ (d x w) 



+ fidy [28yV + {dyW) 



(6) 



(7) 



+ pd x [2d x ud x w + (dyU + d x v)d y w + \Vw\ 2 d x w} 
+ fidy {(dyU + d x v)d x w + 2d y vdyW 

+ \Vw\ 2 d y w} - k(V 2 ) 2 w + fa + V • (PVw), (8) 

where p 2 is the 2D mass density (mass per unit area). The forcing term fa in the right 
hand side (RHS) of (jHJ) tries to bend the membrane upwards or downwards depending 
on the value of the spin at the point (x, y) . a(x,y) is a stochastic variable that flips 
with the transition probability 01]). The term V • (PVw) included in the RHS of (|S]) 
represents the effect of a stress P on the membrane. A homogeneous strain u* x = 7, 
u* y = —X r y/(X + 2/i), u* xy = 0, in a rectangular membrane of length L and width W 
elongated in the x direction and pinned at x = and x = L produces a stress P — E 2 j, 
where E 2 = 4/i(A + p)/{X + 2/i) is the 2D Young's modulus. In addition, the membrane 
may be deflected vertically which yields an additional stress AP = E 2 AS/S, where AS* 
is the change in the membrane area 5* = LW due to bending: 
AS 



S 



^0 



y/l+ \Vw\ 2 - 1 



Then 



4/i(A + fi) 
A + 2/i 



7 + 



LW 



2LW 



dxdy 



2LW 



|Vitf| dxdy. 



w 



Vw\ 2 dxdyj 



(9) 



L L Bonilla and A Carpio 

3. Equilibrium configurations and buckling phase transition 

3.1. Effective free energy provided by spins 

In equilibrium the part of the partition function due to spins is 

Rising = ^2 exp 



f f a(x,y)w(x,y) dxdy 



k B T 



J^exp 
Y[ 2 cosh 



exp 



6 

fa 2 w(x,y) 



<t x,y 



fa 2 o-(x,y)w(x,y) 



exp In < 2 cosh 



x,y 



fa 2 w(x,y) 



In 



2 cosh 



fa?w(x,y) 
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dxdy 



To obtain this equation we have replaced /(•••) dxdy by a 2 Y^ x wherever convenient. 
Added to (JTJ) , the previous equation gives an additional contribution to the free energy 
of the membrane so that its effective free energy is 

/ fa 2 w(x,yY 



"cft : 



— 2„..\2 i ^„,2 i ,,„,2 ^ 



2 cosh 







<ix dy. (10) 



2 (V w) + ^u u + fiu ik - ^ 

The equilibrium configuration is obtained by finding the minimum of F e s with respect 
to the displacement vector. The components u, v satisfy and ([7]) with zero inertia 
(P2 — 0), whereas (|S]) should be replaced by 



A d x 



d x u + d y v + 



\Vw\ 



d x w 



+ Xd v 



d x u + d v v + 



\Vw\ 



OyW 



+ /j,d x \2d x ud x w + (d y u + d x v)d y w + \Vw\ 2 d x w\ 
+ /id y {(d y u + d x v)d x w + 2d y vdyW + \ Vw\ 2 d y w} 

' fa 2 w" 



/t ( V 2 ) 2 ?/; + / tanh 







PV z w = 0, 



(11) 



The flat sheet u = v 
leads to the equation 



w 



is obviously a solution. Linearization of ( fTTT) about it 



-f 2 2 

- k (V 2 ) 2 5w + £-£-6w + P V 2 Sw = 0, 



(12) 



where Pq = E-fy. The flat sheet solution will lose stability at temperatures at which 
Equation (fT2|) has a nonzero solution. Let us assume that the rectangular sheet is simply 
supported at all its boundary points. Then we have Sw = d 2 Sw = at x = and x = L 
and Sw = d 2 Sw — at y — 0, W, so that ^y^"-'" 1 ) = sin(n7rx/L) sm(rrnry/W), and 
therefore [P V 2 - k (V 2 ) 2 ]<^ (n,m) = -P 7T 2 [(n/ L) 2 + (m/W) 2 } - Kit%n/L) 2 + (m/W) 2 ] 2 . 
Substituting in ( fl2l . we find the critical temperatures 

fa 2 



0,, 



^ 2 (xl + wl) 



KIT 



\ L 2 -r w2 ) 



(13) 
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The maximum temperature 6 C = 6^1 gives the critical temperature above which the flat 
graphene sheet is stable. Below this temperature stable ripples form. The same result 
will be obtained below from the macroscopic equations for the average displacement 
vector. 



3.2. Macroscopic equations and critical temperature 



In time-resolved experiments, each framework in a movie consists of data taken during 
one second [30]. Correspondingly, the stochastic equations of motion (jHJ)-(|EJ) with 
the Glauber dynamics (H])-© have to be integrated over long time intervals and 
averaged over shorter intervals during one second. These time averages smooth out the 
displacements of the carbon atoms which become very close to their mean values. The 
equations for the average displacements obtained from (jSJ)-© have to be supplemented 
by the equation for (a(x,y)) = q(x,y) [2Tj . 

' fa 2 w(x,yY 



d t q=- 



^tanh 



e 



Q 



(14) 



To analyze the average equations, we split u = u + Au, etc., in their averages u = (u) 
and fluctuations Au. For a large enough sheet, we ignore the fluctuations and then the 
average values satisfy 



p 2 d t u = X d x I d x u + d y v + 



\Vw\ 



w 



+ /i d y (d y u + d x v + d x wdi 

d 2 t v = A d y (^d x u + d y v + ' 

+ \x d x (d y u + d x v + d x wd y w) , 

\Vw\ 



p 2 d 2 t w = X d x 



d x u + d v v + 



+ pd x [2d x u + (d x w) 



+ jj,d y [2d y v + (d y w) 



(15) 



(16) 



+ Xd v 



d x u + dyV + 



{d x w) 2 + {dyw) 



OyW 



'v- ■ 2 

+ p:d x [2d x ud x w + (d y u + d x v)d y w + iVwl 2 ^-^} 
+ fidy I (d y u + d x v)d x w + 2dyvdyW + \ Vw\ 2 dyw} 
Ap(X + fi) 



K(V 2 ) 2 w + fq + 



A + 2/i 



7 + 



d t q 



a 
2 



tanh 



fa 2 w(x,y) 



- q 



These equations have the trivial solution u = v 
of variables in the linear stability problem (u - 
equations (P = -^27) 



w 



2£wl I l^l 2 «!/)v%,(17) 

(18) 

q = (flat sheet). Separation 



e ip u , ...) leads to the eigenvalue 



p 2 A 2 ij u = (A + 2p)d 2 J) u + (A + p)dl y ip v + pd 2 yip u 



p 2 A 2 ^ = (A + p)d 2 ip u + (A + 2p)d 2 i/j v + pd 2 x i) v , 



(19) 
(20) 
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p 2 A 2 ^ = M + AV 2 ^ - k (V 2 )%,, (21) 

(22) 



2 V 

Eliminating tp q from the last two equations, we obtain 

p 2 A 2 ij w = P V 2 ^ - k (V 2 ) 2 ^ + /* a2 2 f™ . (23) 

As in the previous section, we shall assume that the rectangular sheet is 
simply supported at all its boundary points. We again get eigenf unctions i/jw' mS> = 
sin(mrx/L)sm(rrmy/W), with [P V 2 - « (V 2 ) 2 ]^ n,m) = -P 7T 2 [(n/ Lf + (m/W) 2 } - 
K7r 4 [(n/L) 2 + (m/W) 2 } 2 . Substituting in (123]) . we find the eigenvalue equation 

We have A = at the temperatures n)m given by (1T3"]) . The largest such temperature 
is C = 6>i ;1 : 

= 1 (25) 

° - 2 P„LHl + |i)+^ 2 (l + ^) 2 ' 
Above 9 C , the flat graphene sheet is stable and spins are unpolarized, whereas long non- 
homogeneous ripple patterns with wavelengths related to L and W are formed below 9 C 
and the spins become polarized. Note that increasing the strain 7 decreases 9 C which 
favors ripple destruction. The spins phase transition at 9 C is related to a buckling 
transition in the graphene sheet: The stable state for 9 just below 9 C corresponds to 
a curved sheet all whose atoms are above or below the horizontal plane. As 9 further 
decreases, other critical temperatures # n>m are reached for which patterns corresponding 
to eigenfunctions tfjw' m ^ with n — 1 (resp. m — 1) zeros in the horizontal (resp. vertical) 
direction may appear. These patterns are long ripples in the graphene sheet. For 
a simpler model of a harmonic oscillator coupled to Ising spins undergoing Glauber 
dynamics, the connection between the spins phase transition at 9 C and the instability 
of the rest state of the oscillator and its dynamics are discussed in [311 [32] . The case 
of rippling in a chain of oscillators each connected to an Ising spin undergoing Glauber 
dynamics is considered in |24j. Similar systems with other short-range oscillator-spin 
couplings were studied in relation with the collective Jahn- Teller effect [331 El] • Those 
early studies were focused on whether spin degeneracy is broken below a certain critical 
temperature and on the spin dynamics, but they did not consider the possible patterns 
of the oscillator systems as [23] and the present work do. 



4. Numerical simulations of periodized discrete elasticity. 

We discretize the equations of motion on a hexagonal lattice using the same notation 
as in [261 ESI [27] . Assigning the coordinates (x, y) to the atom A in sublattice 1 as in 
Figure [TJ the three nearest neighbors of A belong to sublattice 2 and their cartesian 
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Figure 1. (Color online) Neighbors of a given atom A in sublattice 1 (dark blue). 



coordinates are ni, ri2 and 713 below. Its six next-nearest neighbors belong to sublattice 
1 and their cartesian coordinates are n^, % = 4, . . . , 9: 

fa a \ fa a 



r» 2V3J V 2 2V3 

f a \ ( a ay/3 

n>3 = I x, y + -j= I , n 4 =\x--,y 

, a ay/3 
n 5 = [x + -,y- — 




n 7 = (x + a, y), n 8 = I x - 

, a ay/3 \ , . 

»■■■= \ x + r y+ ^J- (26) 

In Fig. [U atoms and n 7 are separated from A by the primitive vector ±a and atoms 
n4 and rig are separated from A by the primitive vector ±b. Instead of choosing the 
primitive vector ±b, we could have selected the primitive direction ±c along which 
atoms ns, A and 715 lie. Let us define the following operators acting on functions of the 
coordinates (x, y) of node A: 

Tu = - u{A)\ + [u{n 2 ) - u{A)\ 

+ [u(n 3 ) - u(A)} ~ (flgtt + aj«) j, (27) 

#u = [u(ne) - u(A)] + [u(n 7 ) - u(A)] ~ a 2 <^w, (28) 
Dim = [u(n 4 ) - w(A)] + [u(n g ) - u(A)} (29) 
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D 2 m = [«(n B ) - + [w(n 8 ) - u(A)] 
J^-y + ~ d' z y u I a 2 

A/jM = w(n 7 ) - w(A) ~ (d x u) a, 
A„w = u(n 3 ) — u(A) ~ (c^m) 



a 

7!' 



= [Tw(ni) - Tw(A)} + [Tw(n 2 ) - Tw(A)} 
+ [T«;(n 3 ) - r«;(A)] ~ ^ + ^) V 



(30) 

(31) 
(32) 

(33) 



with similar definitions for points A in sublattice 2. Taylor expansions of these finite 
difference combinations (T, H, D\ and D 2 were defined in [26]) about (x,y) yield 
the partial derivative expressions written above as a — >• 0. We now replace Hu/a 2 , 
(4T — H)u/a 2 and (.Di — /^^/(v^a 2 ), Ahiu/a, \/3A v w/a and 16Sw/a 4 instead of 
<9 2 w, d 2 u, d x d y u, d x w, d y w and (V 2 ) 2 u>, respectively, in (j6]), ([7]) and ([8]) with similar 
substitutions for the derivatives of v and w. We obtain the following equations at each 
point of the lattice: 



p 2 a 2 d 2 u = 4/i Tu + (A + fj,) Hu + (D x - D 2 )v 

v3 

+ ^-t^ [A fe w #w + A v w (A - D 2 )w] + — A^w Tw, 
a a 

p 2 a 2 d 2 v = 4(A + 2/i) + (D 1 - D 2 )u 

v3 

ZL Px 

- (A + ^)/Tv + -^-(A + 2p)A v w Tw 
a 

+ [A ft w (Di - L> 2 )«7 - 3A v w Hw], 



(34) 



2 oz 

p 2 a a t w 



A + 2/i 



(35) 



2A„w A h w 

.HM H (_L>i — lJ 2 )u7 H nw 



+ 



v/3(4T - H)v + ^^(4T - ii> 



a 

A v w 



A h w 



a V3a 
+ ^^[V3(4T - fl> + y/3Hv + {Dx - D 2 )v 

(X 

2A h u 



+ 



a 



(2AT + pH)w + 2v ^ A - v [2\Tw + /i(4T - H)w] 



a 



u9 (a i A h ^ (Di-D2)w , (A h w) 2 + (A v w) 2 /nArr 

+ 2/i A^mH ^ 1 - (2Ai +fiH)w 

V V 3 / a a 2 



+ —Tw(AT - H)w + 4PTw + /a 2 a. 



(36) 
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X X 

Figure 2. (Color online) Ripples formed in a suspended 80 x 80-hexagon graphene 
sheet after (a) 5, (b) 10, (c) 15, and (d) 20 nondimensional time units. Parameter 
values as that in Table[T]but with 5 = ua^fpj \J\ + 2/j = 10~ n . 



Here 

p =^(^^i^r+^n} (37) 

and N = Y2 X y 1 i s ^ ne total number of atoms in the graphene sheet. The spin variable a 
stochastically flips with probability (BJ. Possible defects inserted in the graphene sheet 
are the cores of dislocations. To account for them, we have to write these equations 
of motion in primitive coordinates and periodize all difference operators appearing in 
them along primitive directions. Note that we can use the equations of motion (34\)- $3h}) 
without periodization if there are no defects in the graphene lattice. These equations of 
motion become those in [26] for w = 0. 

fa jT fi 16k f _ aa VP^ 

X+2fi fa 3 \+2fi {X+2n)a 2 ^/X+2~JI 

4.141 x lfr 4 0.4593 0.4428 0.1271 4.084 x 10~ 16 



Table 1. Dimensionless parameters corresponding to f — 3.78 meV/A 3 . 
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Figure 3. (Color online) (a) Ripples formed in a suspended 160 x 40-hexagon 
graphene sheet, (b) Same for a 80 x 80- hexagon sheet. Parameter values are from 
[3"5] at 300K with / = 0.378 meV/A 3 (ten times smaller than that in Table p. 
S = aoy/pjy/x + 2/x = 10~ 4 . 

In our numerical calculations, we have used the values of the 2D Lame moduli at 
300K given in [35J (that agree with experimental measurements p2]), A + 2/i = 22.47 
eV/A 2 , /i = 9.95 eV/A 2 , k = 1.08 eV (see Figure 7 of [36]) and l/a = 27 s (see [26]). 
To calculate /, we balance the energy of a ripple of height wo, E 2 wo\ r (A r is the ripple 
wavelength and E 2 = 4/i(A + /i)/(A + 2ji) the 2D Young modulus), with the coupling 
energy due to spins, fw LW. We find / = E 2 X r /(LW). For a unstrained small sample 
with L = W = 80a and A r = 6.6 A, / = 0.00378 eV/A 3 , the critical temperature is 
11900 times the ambient temperature, and it is 9160 times the ambient temperature if 
there is a 1.5% strain. The critical temperature is much larger for typical samples used 
in the experiments. Thus all graphene samples are in a subcritical buckled state. 

Dimensionless parameters for our model have the numerical values listed in Table 
CD u, v, w are zero for all atoms at the boundaries and for the atoms outside the sheet 
whose displacements enter ([H])-([Hj). In our tests, the average value of \Vw\ is between 
0.1 and 0.001 and it increases with /. The time scale depends on a. To speed up 
simulations, we have used a larger a (and therefore a larger 6 = aa^fp j \J\ + 2/x) than 
in Table [TJ The results are similar if we use smaller a, except that it takes more time 
to arrive at the rippled states shown in the figure. Figure [2] shows the ripples formed in 
a square suspended graphene sheet having 80 x 80 hexagons (12800 atoms) for different 
times. Initially the sheet is planar and the spins are in a random configuration. Figures 
[21 a) and (b) show that there appear small domains in which all atoms are displaced 
vertically upward or downward from the horizontal. These are the ripples. After a 
short transient, they become larger and taller and change very little as time evolves (cf. 
Figures [2(c) and (d)), only atoms in the periphery of a ripple domain flip. Figures [3(a) 
and (b) show ripples in rectangular (160 x 40 hexagons i.e., 12800 atoms) and square 
(80 x 80 hexagons, also 12800 atoms) suspended graphene sheets. The dominant ripple 
wavelength is about 20 lattice spacings (5 nm) in short samples of 80 x 80 spacings 
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Figure 4. (Color online) (a) Ripples formed in a suspended 40 x 40-hexagon square 
graphene sheet containing a vacancy, (b) As in (a) but containing a divacancy. 
Parameters are as in Figure [3] except that 6 = 10 ~ 6 . 



and it goes up to about 45 spacings (11 nm) for longer 160 x 160 samples. There are 
ripples of shorter wavelength superimposed to the longer ones. Over long periods of 
time, significantly larger than the mean time between spin flips, the ripples change. 
This scenario agrees with the observations of Bangert et al [8] . 

Spontaneous ripple generation is caused by coupling with Ising spins. If / = (no 
coupling), the number and appearance of ripples depends on the initial condition. If 
the initial condition compatible with the boundary conditions contains a given number 
of ripples, they persist for subsequent times under the equations of periodized discrete 
elasticity. Adding tension P flattens the ripples somewhat without suppressing them. 
When coupling with Ising spins is restored, ripples such as those in Figures [2] and [3] 
appear after a transient stage even from an initially random state. The average ripple 
height increases with increasing / and decreasing 9. For the boundary conditions we 
use, tension has little effect on the ripples. 

Figure H] shows rippled states for suspended sheets that contain stable defects in 
them: a vacancy in Fig. H^a) and a divacancy or pentagon-octagon-pentagon defect in 
FigfJJb). Both these defects are the cores of dislocation dipoles of equal and opposite 
Burgers vectors and we have calculated them by extending our previous procedure 
from planar graphene [26J to the present model that allows off plane displacements. 
Vacancies and divacancies seem to surf on the ripples affecting little the local curvature 
of the graphene sheet in seeming agreement with experiments [8]. In contrast to this, 
a single pentagon-heptagon effect (dislocation with nonzero Burgers vector) strongly 
affects the sheet local curvature: The sheet is pushed upwards near the pentagon and 
downwards near the heptagon. These effects are due to the coupling with the spins: 
there is no systematic bending of the graphene sheet near defects if / = and ripples 
are introduced through an initial condition. In this case, the initial condition determines 
the local bending of the sheet near defects. 

We have also explored other couplings between spins and carbon atoms. For 
instance we could have coupled the spins to the sheet curvature thereby choosing 
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instead of ([2D- This is similar to the electron-phonon coupling in [TU]. In this case the 
sheet also buckles upward or downward for temperatures below a certain critical value. 
However, simulations of this model show buckling but no ripple formation. 

5. Conclusions 

We have proposed a mechanism to spontaneously produce ripples in a suspended 
graphene sheet. Besides the membrane free energy, carbon atoms are coupled to Ising 
spins at each lattice point that flip randomly according to Glauber dynamics. The 
spins generate forces that pull atoms out from the planar configuration of the graphene 
sheet. After a short transient, there appear domains in which all atoms are vertically 
displaced above or below the horizontal plane. These domains acquire a certain size 
and form ripples with no preferred direction. Changes to the domains are slow and 
are produced only through changes in the atoms at their boundaries. When the spin 
flipping rate is much lower than the frequencies of the membrane, the ripples are long- 
lived metastable states for any temperature. Numerical simulations show that, after 
a transient stage, stable rippling appears even from random initial states. If there 
is no coupling to spins, ripples exist only if they are generated by initial conditions. 
Pentagon-heptagon defects corresponding to dislocations with nonzero Burgers vector 
strongly affect the local curvature of the sheet whereas defects such as vacancies or 
divacancies (corresponding to dislocation dipoles with zero overall Burgers vector) surf 
on the ripples. 
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